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Abstract. Nonlinear tranlational symmetric equilibria with up to quartic flux terms 
in the free functions, reversed magnetic shear and sheared flow are constructed in two 
ways: i) quasianalytically by an ansatz which reduces the pertinent generalized Grad- 
Shafranov equation to a set of ordinary differential equations and algebraic constraints 
which is then solved numerically, and ii) completely numerically by prescribing 
analytically a boundary having an X-point. The equilibrium characteristics are then 
examined by means of the pressure, safety factor, current density and electric field. 
For flows parallel to the magnetic field the stability of the equilibria constructed is 
also examined by applying a sufficient condition. It turns out that the equilibrium 
nonlinearity has a stabilizing impact which is slightly enhanced by the sheared flow. In 
addition, the results indicate that the stability is affected by the up-down asymmetry. 



1. Introduction 

Sheared flows play a role in the transitions to improved confinement regimes in 
magnetic confinement devices, as the L-H transition and the formation of internal 
transport barriers (ITBs), though understanding the physics of these transitions remains 
incomplete. In particular magnetohydrodynamic equilibria with flow, which is the basis 
of stability and transport studies, have been constructed as solutions to generalized 
Grad-Shafranov equations, e.g. Eq. ([I]) below, pQ-pJl]. In connection with the present 
study we refer to our recent contribution [19] in which up-down symmetric nonlinear 
two dimensional cylindrical equilibria with incompressible flow pertinent to the L-H 
transition were obtained. Equilibria relevant to the L-H tranistion usually have peaked 
toroidal current density profiles and safety factors increasing monotonically from the 
magnetic axis to the plasma boundary. 

A necessary requirement for tokamak operation in connection with the ITER and 
DEMO projects is a constant toroidal plasma current, which produces the poloidal 
component of the magnetic field. Among the different options for such non-inductive 
current drive (e.g. electron cyclotron current drive, neutral beam current drive, 
bootstrap current) only the bootstrap current can produce a sufficiently large amount 
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of toroidal current in big tokamaks. The amount of bootstrap current is proportional 
to the pressure gradient. Typically the maximal pressure gradients are located off-axis 
thus leading to hollow current profiles in the plasma associated with reversed magnetic 
shear. Static equlibria with reversed magnetic shear was the subject of [20j-[26j. 

The stability of fluids and plasmas in the presence of equilibrium flows non parallel 
to the magnetic field remains a tough problem reflecting to the lack of necessary and 
sufficient conditions. Only for parallel flows few sufficient conditions for linear stability 
are available [21]- [29]. In previous studies we found that the stability condition of |29j 
is not satisfied for the linear equilibria of [10] and [15] while it is satisfied within an 
appreciable part of the plasma for the nonlinear equilibria of p2], [18] and [19]. This 
led us to the conjecture that the equilibrium nonlinearity may act synergetically with 
the sheared flow to stabilize the plasma. 

Aim of the present study is to extend our previous paper [Tjj] in two respects: up- 
down asymmetry and reversed magnetic shear. As in [19] non magnetic field aligned 
equilibrium flows will be included. In this respect it is noted that a synergism of reversed 
magnetic shear and sheared poloidal and toroidal rotation, consisting in that on the 
one hand the reversed magnetic shear plays a role in triggering the ITBs development 
while on the other hand the sheared rotation has an impact on the subsequent growth 
and allows the formation of strong ITBs, was observed in JET [30] and DIII-D [3 1 J . 
In addition here the above conjecture about a combination of stabilizing effects of 
equilibrium nonlinearity and plasma flow will be checked. The reason for considering 
translational symmetry is the many free physical and geometrical parameters involved 
in connection with the flow amplitude, direction and shear, equilibrium nonlinearity, 
symmetry and toroidicity. Thus, in the presence of nonlinearity one first could exclude 
toroidicity. 

The organization of the paper is as follows: In the second section we briefly 
present the general setting for the translatinally symmetric equilibrium equations 
with incompressible flow and introduce the ansatz reducing the problem to a set of 
ordinary differential equations (ODEs) and algebraic constraints. In section 3 up-down 
asymmetric equilibria are constructed quasianalytically and their characteristics are 
studied. In section 4 equilibria with a lower X-point are derived numerically by imposing 
analytically the boundary shape. A stability consideration of the equilibria obtained is 
made in section 5. Section 6 summarizes the conclusions. 

2. Translational symmetric equilibria with flow 

The equilibrium of a cylindrical plasma with incompressible flow and arbitrary cross- 
sectional shape satisfies the generalized Grad-Shafranov equation [I], [7] 



for the poloidal magnetic flux function if). Here, M p (if)), P s (4>), p(if>) and B z {ifj) are 
respectively the poloidal Alfven Mach function, pressure in the absence of flow, density 




(1) 
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and magnetic field parallel to the symmetry axis z, which are surface quantities. Because 
of the symmetry, the equilibrium quantities are ^-independent and the axial velocity v z 
does not appear explicitly in Eq. (JT]). SI units are employed unless otherwise stated 
(see section 5). Derivation of Eq. ([1]) is based on the following two steps: First express 
the divergence free fields in terms of scalar quantities as 

B = B z Vz + Vz xVip 
m = V 2 i)Vz -Vzx VB Z 

pv = pv z Vz + Vz x VF (2) 

The velocity v relates to the electric field, E = — V<3? (where is the electrostatic 

potential), by Ohm's law, E + v x B = 0. Second, project the momentum equation, 
p(v • V)v = j x B — VP, and Ohm's law, along the symmetry direction z, B and Vip- 
The projections yield four first integrals in the form of surface quantities (two out of 
which are F(tp) and <£>(■?/>)), Eq. (CQ) and the Bernoulli relation for the pressure 

P = PsW - ^M p 2 (V0lW| 2 (3) 

Because of the flow P is not a surface quantity. Also the density becomes surface 
quantity because of incompressibility and M 2 (-?/>) = (F' (ip)) 2 /(pop). Five of the surface 
quantities, chosen here to be P s , p, B z ,Mp and v z , remain arbitrary. 
Using the mapping 

r#) = f[l - M 2 p (g)Y/ 2 dg, (M 2 < 1) (4) 

•/ 

Eq. ([1]) is transformed to 

= (5) 




du 



Note that transformation (jlj) does not affect the magnetic surfaces, it just relabels them. 
Eq. (jSj) is identical in form with the static equilibrium equation. 

In the present study we assign the free function term in Eq. (jSJ) as 

( B 2 \ u 2 u 3 u A 

+ -y j = c + c t u + c 2 — + C 3 y + C 4 y (6) 

to obtain 

u xx + u yy + Ci + c 2 u + c 3 u 2 + c^u 3 = (7) 

where (x, y) are the usual cartesian coordinates. The form of this equation leads us to 
introduce the following up-down asymmetric ansatz for the flux function which enables 
reduction of the equilibrium problem to a set of ordinary differential equations and 
first-order constraints: 

_ N 2 {x)y 2 + N 1 {x)y + f{x)D {x) 
U y 2 + D 1 (x)y + D (x) U 
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This is an extension of the respective ansatz for up-down symmetric equilibria we 
introduced for the first time in [19]. Inserting Eq. (JSJ) into (J7J), after a rather lengthy 
calculation the latter is transformed into a fraction (F), the nominator of which is a 
polynomial of y of sixth order. Equating this nominator to zero, from the y 6 -term we 
obtain 

N'z + ci + c 2 N 2 + c 3 N 2 + c 4 iV 2 3 = (9) 

From the 7/°-term it follows 

2(N 2 - f) _ -WxjNi - /A 

From the y 4 , y 5 -terms one yields 



/ + " n JJ V n2 J ' +c 1 + c 2 f + c 3 f 2 + c 4 f i = (10) 



< = [-iv 2 (iv 2 - /)d + (/D 2 (iV 2 - /)- 

-D (JVi - /^i)(iVi + N 2 Dx))G 5 + £> iV 2 (^i - /A)C 4 ] (11) 

£>i = ^ [-(iV 2 - f)G 1 + (D 2 (N 2 -f)- 2D Q D 1 (N 1 - fD 1 ))G 5 

+D (N 1 - fD 1 )G 4 \ (12) 



' \{Nx - N 2 D l )Gi - DKN, + N 2 D X - 2fD x )G 



DD 

+D%(N 2 - f)G 4 ] (13) 

where DD = D (Nx - /£>i)(JVi - N 2 D X ) + D 2 (N 2 - f) 2 and the functions G u G 4 , G 5 
are given in the Appendix. From the y 3 , |/ 2 -terms we obtain respectively the first-order 
constraints Cx, C 2 

(D? + 2D O )<(0>') - (N 1 + N 2 D l -2fD l )D'^<f ) ,<p') - 

- (N 2 D + NxDx + fD )Dl((/>, </>') + G 3 = (14) 

and 

2D DxN';(<f>,<f>') - (D (N 2 -f) + Dx(Nx - fD^D'^J) - 

- D (Nx + 0') + G 2 = (15) 

Here denotes collectively all the functions appearing in Eq. (JED, and Eqs. (JUKES]) are 
used. The functions G 2 , G% are also given in the Appendix. Note that Eqs. (TT41 . (floT) 
are not second order "evolution" differential equations but rather first-order constraints 
to be fulfilled during numerical integration. This justifies the introduction of ansatz 
(jHJ) which results in a simple numerical treatment of the equilibrium through ordinary 
differential equations and algebraic constraints in contrast and alternative to the full 
numerical treatment of section 4. 



3. Class of quasianalytic solutions 



The system of Eqs. (!9l fT3]) is integrated numerically using high-precision numerical 
integration with very small step size, due to its extreme complexity and nonlinearity. 
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The magnetic axis regularized with respect to the geometric center is taken to be 
"Shafranov shifted''^ at x a = 1 + x s = 1.1. We take the following ITER-pertinent 
geometrical data: a = 2 m, Rq = 6.2 m for the minor and major radius of the "torus" 
respectively and the inverse aspect ratio is eo = 0.32. The bounds for the x— variable 
are x m i n = 1 — eo and x max = 1 + e . The integration begins from x a forward up to 
x max and backwards up to x min . For the following values of the parameters C\ = 52.0, 
C2 = —0.4, C3 = 0.1 and C4 = 0.1 and for initial conditions N 2 = —1.4, N' 2 = 0.95, 
iVi = -0.15, N[ = -1.05, D 1 = 0.15, D[ = -0.137, D = 1.05, D' Q = 0.1, / = 7.5 
and / = 0.0 we obtained the solution of Fig. 1. The first-order constraints of Eqs. 
( TTIlfToT) were monitored during the integration and they were kept to very low values. 
The product of the constraint C\ with the average value of y 3 (taken to be equal to 0.2) 
and the constraint C 2 with the average value of y 2 , that appear in the nominator of the 
final fraction (F) (see the text just below Eq. (|H])) were kept bounded to \C\\ < 0.07 
and [C2I < 0.18. This combined with the fact that the denominator of this fraction 
has a positive definite value greater or equal to 1.0, is an additional argument that this 
fraction is very close to zero and that the presented equilibrium is indeed an acceptable 
solution of Eq. (J7j). 

The bounding surface (shown in green) corresponds to Uf, = 5.5 Wb while the 
magnetic axis to the value u a = 7.8 Wb. The magnetic axis is located at (x a ,y a ) = 
(1.1, —0.07). A quartic fitting yields expressions for the functions N 2 , N\,Di, D , f such 
as the following equation 

N 2 = 5.527a; 4 - 23.78a; 3 + 11.81a; 2 + 31.95a; - 27.27 (16) 

We stress however the fact that for a correct representation and plotting of the various 
functions all the higher order expansions and more precise numerical parametric values 
are needed. Plotting Eq. ([8]) using Eq. ( {TBI) will not yield the correct result of Fig. 1, 
which occurs through the precise numerical results for these functions. 
The MHD safety factor is defined as [32] 

_ djltor/dV 1 djltor/dV . . 

q[U) ~ di) pol /dV ~ 2tt dij/dV [ ' 

Using dV = 2ir\J\dil ! d8, J -1 = W • (V</> x V^) and expressing the fluxes in terms of the 
magnetic field, one can cast ([17]) in the form of a line integral on each constant— u curve. 
The detailed evaluation of q for the present equilibrium of Fig. 1 yields the curve of 
Fig. 2 with strong reversed magnetic shear. For this result it was used cq = 200.0 
in Eq. (J6]). Also for the axial magnetic field it was adopted the typical tokamak 
diamagnetic function B z = B z0 (l + 7(1 — ^)), shown in Fig. 3, with 7 = 0.1 and 
B z0 = 3.2 T . Then, the static pressure function P s (u) is computed by Eq. ([6]), while 
for the flow function in Eq. (0} it was used M 2 = M pa (^ - l) 2 - 5 with M pa = 0.1 and 
Ub = 5.5 Wb < u < u a = 7.5 Wb. The pressure (Eq. ([3])) is shown in Fig. 4, normalized 
to its center value of Po = 2.1046 x 10 5 Pa. Also instead of the axial velocity v z , the 

f The term "Shafranov shift" here means that the equilibrium in addition to up-down is left-right 
asymmetric and is not connected to the toroidicity which vanishes in cylindrical geometry. 



Two dimensional nonlinear cylindrical equilibria. 



6 



corresponding Mach function M| is chosen similar to the poloidal one (M| ~ Mi 2 ) with 
M za = l.lM pa . 

The electric field for equilibrium of Fig. 1 is shown in Fig. 5. Here the choice 
P = P a (it ~~ ^ as ^° een ma de for the density with p a = 4.0 x 10~ 7 Kgr m~ 3 . The 
maximum of E increases with the flow parameter M pa but the position of the maximum 
is not affected by the flow in agreement with the results of [TJ [16]. The hollow axial 
current density profile in the midplane y — is shown in Fig. 6 in consistence with the 
negative magnetic shear curve for the safety factor of Fig. 2. 




X 



Figure 1. Equilibrium for the initial conditions given in the text of section 2. The 
bounding surface, shown in green, corresponds to ui, = 5.5 Wb while at the magnetic 
axis, which is located at (x a ,y a ) = (1-1, —0.07), u a = 7.8 Wb. 



4. Numerical asymmetric equilibrium with X-point 

We consider now the direct numerical solution of Eq. ((7j) with a prescribed boundary 
possessing a divertor null X-point. We first will specify the boundary. The boundary 
conditions for the flux function u is Ub = 1 Wb on the prescribed boundary curve and 
u a = Wb on the magnetic axis. The magnetic axis is taken to be Shafranov shifted at 
x a = 1 + x s = 1.1. The model thus has three free parameters which are the Shafranov 
shift x s , the elongation n and the triangularity 5. Their values are taken as x s = 0.1, 
k = 1.86 and 5 = 0.5, in accordance with the corresponding data of the ITER project. 
The bounding flux surface, for the aymmetric case, is shown in Fig. 7. We take the 
values a = 2 m, Rq = 6.2 m for the minor and major radius of the "torus" respectively 
for which the inverse aspect ratio is eo = 0.32. 

The equation for the upper part of the bounding flux surface, which if taken to 
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5.5 6 6.5 7 7.5 



u 

Figure 2. The safety factor for the equilibrium of Fig. 1 presenting a strong negative 
magnetic shear region. The outer bounding surface corresponds to u — Ub = 5.5 Wb 
while at the magnetic axis u = u a = 7.5 Wb 

1 .005 1 1 1 1 1 1 1 




x 



Figure 3. The axial magnetic field B Zl for equilibrium of Fig. 1 normalized to its 
center value of B z q — 3.2 T. 

hold for the lower part as well would give a symmetric bounding surface, is 
Xb = 1 + €qcos(t + asin{r)) 

Vb = y m axsin{r) (18) 

where y max = ^eo with 5 = (1 — xs)/e , and a = sm _1 (J). Thus the following relations 
hold: xs — 1 — Seo and 0s = n — tan~ 1 (n/5). The parameter r is any increasing function 
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Figure 4. The pressure for equilibrium of Fig. 1 normalized to its center value of 
P = 2.1046 x 10 5 Pa. 




1.2 1.3 



Figure 5. The electric field E z , for equilibrium of Fig. 1 normalized to its maximum 
value of E = 5.4851 kV m _1 . 



of the polar angle 9, satisfying r(0) = 0, t(tt) = n and r(8s) = ir/2. In our model we 
take 

t{9) = t 9 2 + h6 n 

° 7T0? " ^TT"- 1 
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Figure 6. The axial current density at the midpalne A(x, y = 0) for the equilibrium 
of Fig. 1. It is hollow in connection with the negative magnetic shear of Fig. 2. 




Figure 7. Bounding flux surface for the asymmetric case defined by Eqs. (|18ti21[) . 
possesing a divertor, null X-point at Px(xx,Vx) = (0.9139, —0.6105). 

-9 S + 2 fl9l 

1 7T0? - 

with n = 8. In order to complete the asymmetric bounding curve we specify now the 
lower part of it (y < 0) as follows. The left lower branch of the curve is given by 

xi, = 1 + 60005(6*) 

y b = -[2 Pl e Q (l + cos6)] l l 2 



Two dimensional nonlinear cylindrical equilibria. 



10 




Figure 9. The stability function A for the equilibrium of Fig. 1. At the upper part 
(y > 0) of the equilibrium it mostly assumes positive values, while at the lower part 
(y < 0) it assumes negative values. 



= o (fr ay (?r - 9 - 2n - e&) (20) 

2e (l + cosOs) 



while the right lower branch of the curve is given by 

Xb = 1 + €ocos(9) 

y b = - [2p 2 e (l - cosO)] 1 / 2 
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7.5 - 



7 - 



6.5 - 



5.5 
-0.8 




0.6 



Figure 10. The poloidal magnetic flux function u(x — 1.1, y) for the equilibrium of 
Fig. 1. It is slightly up-down asymmetric with repsect to y = 0. 



1 .5 r 



c4=0.10 
c4=0.15 
c4=0.20 



0.5 - 



-0.5 - 



0.9 



1.1 

x 



1.2 



1.3 



Figure 11. The stability function at the midpalne A(x,y = 0) for the equilibrium of 
Fig. 1, for various values of the nonlinear constant C4. It appears that nonlincarity 
acts in favour of the stability. 



P2 



v 2 

Umax 



(2tt - 9 S < 9 < 2tt) 



(21) 



2eo(l — cos9$) 

The divertor null X-point is located at xx = 1 + e cos9s = 0.9139 and yx = —y ma x = 
-0.6105. 

The Laplacian operator u xx + u yy is discretized on a rectangular grid where we have 
(1 — e ) < x < (1 + eo) and —y rna x < y < ymax with grid step h. The nine point formula 
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Figure 12. The stability function at the midplane A(x,y = 0) for the equilibrium 
of Fig. 1, for various values of the flow parameter M pa . It appears that flow acts in 
favour of the stability, though in a weak manner. 

is then employed [33] 

^ 2ui ' j = Qhp \ Ui + 1 d+ 1 + 4M *+bi + u i+i,j-i + 4u id+1 + 4Mij_i + 

+4u^ hj + iti-u-i - 20u id ] (22) 

and when substituted into Eq. ([7j), the latter is written as uff = G(u)\°j d ^ and is 
solved iteratively. The last term of Eq. ( )22|) is taken to be the u^ w ^ and the iterations 
stop (i.e. reaching convergence to the solution) when \u\™ ew ^ — uf^ \ < 0.001. We stress 
again that the conditions Uf, = 1 on the boundary and u a = on axis are imposed in 
every iteration. For h = 0.02, and for the following values of the constants of Eq. ([7]), 
Ci = —10.0, c 2 = 2.0, c 3 = 1.1, c 4 = 1.1, a number of N = 165 iterations were needed 
to obtain the desired accuracy. The solution is shown in Fig. 8. Although it seems that 
the solution is dependent on the specific value of the h chosen, from the discrete two- 
dimensional matrix of u it j produced, a two-dimensional Lagrange fitting is performed 
that yields a polynomial in (x,y). This is h- independent. It is noted that the ripples of 
the flux function, appearing near the boundary of the equilibrium are due to numerical 
instabilities that are inevitable present in the calculation. 

5. A stability consideration 

We now address the important issue of the stability of the solutions constructed with 
respect to small linear MHD perturbations by means of the sufficient condition of [29|. 
This condition concerning internal modes states that a general steady state of a plasma 
of constant density and incompressible flow parallel to B is linearly stable to small 
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three-dimensional perturbations if the flow is sub-Alfvenic (M 2 < 1) and A > 0, where 
A is given below by f l23|) . Consequently, using henceforth dimensionless quantities we 
set p = 1. Also, for parallel flows (v = MB) it holds M p = M z = M. In fact if 
the density is uniform at equilibrium it remains so at the perturbed state because of 
incompressibility. In the u-space for axisymmetric equilibria A assumes the form 

A = - g 2 (jxVu)-(B- V)Vw+ 




with 



//■ „ 1 , ^ (23) 
^;h-(m p 2 )'5 2 /2 



M 2 
p 



This condition, although complicated is accurate (a proof is provided in [29]) and all 
the computations and conclusions have been performed with great care. Specifically, its 
application to the equilibria constructed in sections 3 and 4 led to the following results: 

(i) Even a weak up-down asymmetry affects stability as indicated in Fig. 9 where we 
have checked thoroughly the values of the function A for the equilibrium of Fig. 
1. It turns out that at most of the upper part of the equilibrium, where y > 0, A 
assumes small positive values, while for the lower part of the equilibrium, (y < 0), 
it assumes small negative values. This slight A-up-down asymmetry is connected to 
the respective slight asymmetry of the flux function u with respect to the vertical 
position y; the latter can be seen in Fig. 10 where the profile u(x,y) is given at 
the point x — 1 + x s — 1.1. Though A < does not necessarily imply an unstable 
equilibrium because the condition is sufficient, the above result is consistent with 
the fact that up-down asymmetry may make the plasma unstable. For external 
modes this might relate to the vertical instability, e.g. |3"4"] . 

(ii) The non linearity favours the stability as it is shown in Fig. 11 where A is plotted 
as a function of x at the midplane y = for increasing values of the non-linearity 
constant c 4 . 

(iii) The flow has a slight stabilizing effect. An example is given in Fig. 12 in connection 
with the flow parameter M pa . 



6. Summary 

We have constructed and studied nonlinear translational symmetric equilibria with 
strong reversed magnetic shear and sheared incompressible flow non parallel to the 
magnetic field on the basis of a generalized Grad-Shafranov equation (Eq. (DO)). This 
equation can be transformed to one identical in form with the static Grad-Shafranov 
equation which we have solved in a couple of alternative ways: i) by using an ansatz 
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for the unknown magnetic flux function (Eq. (jSJ) which reduces the original equation 
to a set of ODEs and algebraic constraints (Eqs. (I10H15I) ); then this set of equations 
is solved numerically, and ii) fully numerically by prescribing anallyticaly a diverted 
boundary (Eqs. f)18H2ip ). The equilibria constructed are typically dimagnetic (Fig. 3), 
have peacked pressure profiles (Fig. 4), hollow toroidal current densities (Fig. 6) and 
electric fields possessing a maximum (Fig. 5). The maximum of E takes larger values 
as the flow amplitude increases but its position is insensitive to the flow. 

For parallel flows application of a condition for linear stability implies that the 
equilibrium nonlinearity has a stabilizing effect together with a weaker stabilizing impact 
of the sheared flow in agreement with past nonlinear equilibrium studies p21 UHl H~9] . 
Also even a small up-down asymmetry influences stability. 

Finally it would be interesting to try constructing equilibria with flow and reversed 
current density in connection with non nested magnetic surfaces thus generalizing the 
static ones of [2TJ[22j[2l] and extend the study to axially symmetric equilibria by possibly 
generalizing the ansatz (jBJ) in order to examining the impact of toroidicity. 

Appendix: Quantities appearing in the ODEs and algebraic constraints 

pins) 

The functions Gi, G 4 , G 5 appearing into Eqs. (TTT1) . (|T2|) . ({TBI are given by 

G x = - ADoD^Ni -f) + 4D?(JVi - /£>i) - 6D (A 1 — fD — 1) 

- 2D 0j D;(A; - fD[ - f'D l ) + 2(D' ) 2 (N l - fD l ) - 2D 2 D'j' + 
+ Cl DlD x + c 2 DlN l - c-sfD 2 ^ + 2c 3 fD 2 Nl + 

+ 3cJ 2 D 2 N l - 2c 4 fD 2 D l (24) 

G 4 = - 2(A 2 - /) + 2Dl(Nl ~ fDl) - 2D' (N' 2 - /) - 

- 2D;(A; + D X N' 2 - A 2j D;) + 2 Cl D\ + 2c 2 N x D x + c 3 N 2 + 
+ 2c 3 fD N 2 + 2c 3 N 1 N 2 D 1 - c 3 f 2 D - c 3 N 2 Dj - c 3 D N 2 - 

- c 4 fD - caNId{ + 3cJD N 2 + 3c 4 N 2 N 2 - 2c 4 A 2 3 D (25) 

G 5 = - 2D[N' 2 + ciDi + c 2 N, + 2c 3 N 1 N 2 - c 3 N 2 D 1 - 2c 4 N'^D 1 + 

+ 3c 4 A!A 2 2 (26) 

The functions G 2 , G 3 appearing into Eqs. (11411151) are given by 

G 2 =- 2D\{N 2 -f)+ ggl^LZ f D ^ + 2ci D D 2 1 - c 3 fD,D\ - 

- ctfD Dl - 10D (N 2 -f)+ 4D 1 (A 1 - fD x ) - c 3 f 2 D 2 - 

- 2cJ 3 D 2 + 2f'D' Q {Dl + 2Do) - 2D' (f ' D Q + fD' ) 

- 2D 1 D' 1 (f'D + fD' ) - c 3 N 2 D 2 - c 4 N*D 2 - 2N' 1 D' D 1 - 

- 2A;D d D; - 2N' 2 D D' + 2fD (D' 1 ) 2 + 2N 2 {D' ) 2 + AN^D^ + 
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+ 2c 2 N 1 D D 1 + 2c 3 fD D 1 N 1 + c 3 N 2 D + 2c 3 fD 2 N 2 + 3cJD N 2 + 
+ 3c 4 f 2 D 2 N 2 (27) 



and 



G 3 = - AD X {N 2 - /) + JUX) - 2c 3 f 2 D D 1 - 2c 4 / 3 A ) A + 




o 



+ ADJ'D', - 2c 3 N 2 D D 1 - 2c,N 2 i D D 1 - 2(f D + fD' )D[ - 
- 2N[D' - 2N' 1 D 1 D' 1 - 2N , 2 D' D l - 2N 2 D D' 1 + 2iVi(L>i) 2 + 
+ AN 2 D' D[ + 2iVi - 2jV 2 £>i + c x D\ + 201,00^! + c^Dj + 
+ 202^1)0 + ^c 3 fD Q N x + cs^A + 2c 3 fD D 1 N 2 + 
+ + 2c 3 N 1 N 2 D + c A Nf + 6c 4 Ar 1 7V 2 /A) 



(28) 
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